


function mininfo...
         =min_sa(fun, p, theta, v, T, epsilon, N_epsilon, N_S, c, N_T, r_T, l_bd, u_bd);

%**************************************************************************************************************
%
% This routine implements the simulated annealing algorithm to MINIMIZE the function fun(theta).
% 
% The algorithm and most of the notation comes from Corana et al. (1987): "Minimizing Multimodal 
% Functions of Continuous Variables with the Simulated Annealing Algorithm". Transactions on Mathematical 
% Software, Vol. 13, No. 3, 262-280. This article is an excellent description of the algorithm.
%
% Goffe et al. (1994) (J of Econometrics 60) give good insights too, but they leave out a couple of important
% details necessary for the implementation of the algorithm. 
%
% The output of the procedure is a cell array "mininfo". The first cell mininfo{1} contains the minimizer vector 
% and the second cell mininfo{2} contains the minimum value. 
%
%**************************************************************************************************************

% Argument list:

% fun                          The name of objective function to be minimized. Note that the first argument 
%                              of the procedure "min_sa" is a function: "fun" is the name of the M-file in which
%                              the function is defined (this file must be placed in the work directory). When the procedure 
%                              min_sa is called, the input argument "fun" must be preceded by a @.                      
%                                       
% p=2;                         "fun" is a function of p variables 
% theta=[0.1 0.1];             a 1 x p vector, the starting point for the optimization procedure. In general,
%                              the argument of "fun" is denoted by "theta". 
% v=[.5 .5];                   a 1 x p vector, the starting step size
% T=1000;                      the starting temperature
% epsilon=10^(-6);             terminating criterion 
% N_epsilon=4;                 the number of iterations over which the termination criterion 
%                              must be satisfied for code to stop
% N_S=20;                      how many search cycles over each coordinate before adjusting the step size
% c=[2 2];                     a 1 x p vector, a controls step size adjustment for each coordinate
% N_T=100;                     how many step adjustments before reducing the temperature
% r_T=0.85;                    temperature reduction coefficient
 
% l_bd=[-1 -1 ];               a 1 x p vector, the ith component of l_bd is a lower bound for the corresponding 
%                              component of theta (i.e. we search for the optimal theta in a box); could put 
%                              "-inf" if desired.
% u_bd=[ 1  1 ];               a 1 x p vector, each component is an upper bound for the corresponding component 
%                              of theta (i.e. we search for the optimal theta in a box); could put 
%                              "inf" if desired.

i=0;                                 %index of succesive minimizers (not actually used for anything)
k=0;                                 %index of successive temperature reductions (not actually used for anything)

f=feval(fun, theta);                 %This line calculates f(theta). The reason why the "feval" command must be used
                                     %is because "fun" is just the placeholder for an actual m-file.        

theta_opt=theta;                     % initial optimizer=starting point
f_opt=f;                             % initial optimium=initial function value
n=zeros(1,p);                        % the ith component of the 1 x p vector n gives the number of accepted steps in the
                                     % ith coordinate direction over the current search (N_S) cycle

f_star=f*ones(1, N_epsilon+1)+1;     % the vector f_star contains the current and last N_epsilon optimal
                                     % function values (used to test for stopping). Initially, the value of f_star is set 
                                     % so that the stopping condition is not satisfied.
                                     
e=eye(p);                            %p dimensional unit matrix                                     
                                   
while ~ (     (abs(f_star-f) < epsilon)   &   ((f-f_opt) < epsilon)   )
%if terminating condition NOT satisfied, proceed    

    theta=theta_opt; %start at the current optimizer
    f=f_opt; %current optimal function value
    
    disp(['T=' num2str(T) '   theta_min=(' num2str(theta_opt) ')'...
        '   f_min=' num2str(f_opt) '  v=(' num2str(v) ')']) %display some stuff
   

        for t=1:N_T %temperature cycle                                 
    
                                 
            for j=1:N_S %search cycle; N_S=number of searches through coordinates before step size adjustment
        
            
                for h=1:p %examine each coordinate direction; h=current search direction
    
                        theta_new=theta+(2*rand-1)*v(h)*e(h,:);  %new candidate for theta; rand=unif[0,1]
                                                                 %step a random fraction of the step size v(h) in direction h
                                                
                        if ( theta_new(h)>=l_bd(h) & theta_new(h)<=u_bd(h) ) %check if candidate theta is 
                                                                             %in domain (i.e. "search box")
                            
                            f_new=feval(fun, theta_new); %fn value at candidate point
                            
                        else
                            
                            f_new=inf; %if fn value outside domain, set function value=inf-->new point will be rejected 
                            
                        end

       
                            if f_new < f %accept new point if new function value is lower than current value
            
                                    theta=theta_new; 
                                    f=f_new;
                                    i=i+1;       %index of accepted points
                                    n(h)=n(h)+1; %number of acceptances in direction h in the current search cycle
            
                                        if  f_new < f_opt %if new function value < current optimum,
                                                          %redefine optimum
                    
                                                theta_opt=theta_new;
                                                f_opt=f_new;
                    
                                        end
                 
                            else %if f_new >= f, apply Metropolis criterion to decide whether to accept
                                
                                        q=exp((f-f_new)/T); %acceptance probability (gets smaller as T decreases; also,
                                                            %for a given T, smaller if function would increase too much)
            
                                        if rand < q %accept (otherwise reject and start new direction)
                        
                                                theta=theta_new;
                                                f=f_new;
                                                i=i+1;   %index of accepted points
                                                n(h)=n(h)+1; %number of acceptances in direction h
                                                             %in the coordinate search cycle
                        
                                        end
                    
                 
                            end %end of acceptance tests
        
    
                end %h (cycle of coordinate directions): new direction
   
    
        end %j (N_S searches over each coordinate direction before step adjustment)

        % Update the step vector v so that approx 50% of new points are accepted (see Corana et al. 1994)

        for h=1:p
    
                if n(h)>0.6*N_S
        
                      v(h)=v(h)*(1 + c(h)*(n(h)/N_S - 0.6)/0.4 );
        
                end

                if n(h)<0.4*N_S
        
                      v(h)=v(h)*(1 + c(h)*(0.4 - n(h)/N_S)/0.4 )^(-1);
        
                end
    
                %if neither condition is satisfied, the corresponding coordinate remains unchanged
    
         end %h
    
         n=0*n; %reset acceptance counter to zero 
            
    end %t (N_T cycles with the same temperature)            
                
    %Reduce temperature

    T=r_T*T;
    k=k+1; %number of times temperature is reduced; 

    for u=N_epsilon+1:-1:2
    
        f_star(u)=f_star(u-1);
    
    end

    f_star(1)=f;


end %the while cycle with the stopping criterion 

mininfo={theta_opt f_opt};  


